*************************************************************
* The following code is part of the replication archive of  * 
* Genovese, Kern & Martin (2016) "Policy Alternation", ISQ  * 
*************************************************************

cd "/Users/.../gkm_replication_isq/carbon policy"

* This Stata MSTAR code is presented in Franzese, Hays and Kachi (2009). We have normalized co2tax_value.

************ 
* Figure 4 *
************
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear
cd "/Users/.../gkm_replication_isq/carbon policy/Figures"
collapse (sum)  co2tax_impl co2trade_reg_impl, by(year)
twoway  (line co2trade_reg_impl year) (line co2tax_impl year), ylab(0(5)20, nogrid) xtitle("Year") ytitle("Number of countries") graphregion(fcolor(white)) plotregion(fcolor(white))  scheme(s2mono)  legend(label(1 "CO2 Trade") label(2 "CO2 Tax")) saving(Figure4, replace) 
clear

************************** 
* Summary stats: Table 2 *
**************************

use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear
sum co2tax_value_pc co2taxnormal co2tax_value gdppercapita gdppercapitasq co2ktpercapita co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration dpi ratio_over_eutob allowances_gdp_lag
sutex co2tax_value_pc co2taxnormal allowances_gdp gdppercapita co2ktpercapita euintegration dpi ratio_over_eutob energydependencektofoileq govteffectiveness, minmax

********************************************************************************* 
cd "/Users/.../gkm_replication_isq/carbon policy"

********************
* Results: Table 3 * 
********************

/*
* 0. Focal policy (Co2 tax) weighted by geographic distance - no allowances

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1  sigma
tempvar A rSL1  
gen double `rSL1'=`rho1'*SL1
*gen double `rSL2'=`rho2'*SL2
scalar p1 = `rho1'
*scalar p2 = `rho2'
matrix p1W1 = p1*W1
*matrix p2W2 = p2*W2
matrix IpW = I_n - p1W1
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`mu', 0, `sigma'))
end

**********************
*Open Data For Weights
**********************

* pick the following spatial weights
spatwmat using "/Users/.../gkm_replication_isq/carbon policy/distance_gleditch_ind.dta", name(W1)  standardize

**************************
*Open Data for Regression
**************************
drop _all
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear

global nobs = 187  
matrix I_n = I($nobs)
global Y co2taxnormal
global X lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration dpi ratio_over_eutob 

mkmat $Y, matrix(Y)
matrix SL1 = W1*Y
svmat SL1, n(SL1)
*matrix SL2 = W2*Y
*svmat SL2, n(SL2)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X) (rho1:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
*ml init rho2:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo M1
drop SL11
*/

* 1. Focal policy (Co2 Tax) weighted by geographic distance + CO2 allowances

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1  sigma
tempvar A rSL1  
gen double `rSL1'=`rho1'*SL1
*gen double `rSL2'=`rho2'*SL2
scalar p1 = `rho1'
*scalar p2 = `rho2'
matrix p1W1 = p1*W1
*matrix p2W2 = p2*W2
matrix IpW = I_n - p1W1
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`mu', 0, `sigma'))
end

**********************
*Open Data For Weights
**********************

* pick the following spatial weights
spatwmat using "/Users/.../gkm_replication_isq/carbon policy/distance_gleditch_ind.dta", name(W1)  standardize

**************************
*Open Data for Regression
**************************
drop _all
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear


global nobs = 187  
matrix I_n = I($nobs)
global Y co2taxnormal
global X lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration dpi ratio_over_eutob allowances_gdp_lag

mkmat $Y, matrix(Y)
matrix SL1 = W1*Y
svmat SL1, n(SL1)
*matrix SL2 = W2*Y
*svmat SL2, n(SL2)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X) (rho1:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
*ml init rho2:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo M2
drop SL11


/*
* 1.b Alternative policy (Co2 Allowances) weighted by geographical distance w/o focal policy (Co2 tax) weighted by geographical distance

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1  sigma
tempvar A rSL1  
gen double `rSL1'=`rho1'*SL1
*gen double `rSL2'=`rho2'*SL2
scalar p1 = `rho1'
*scalar p2 = `rho2'
matrix p1W1 = p1*W1
*matrix p2W2 = p2*W2
matrix IpW = I_n - p1W1
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`mu', 0, `sigma'))
end

**********************
*Open Data For Weights
**********************

* pick  the following spatial weights
spatwmat using "/Users/.../gkm_replication_isq/carbon policy/distance_gleditch_ind.dta", name(W1)  standardize

**************************
*Open Data for Regression
**************************
drop _all
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear

global nobs = 187  
matrix I_n = I($nobs)
global Y co2taxnormal
global X lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration dpi ratio_over_eutob allowances_gdp_lag
global Z allowances_gdp_lag

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL1 = W1*Z
svmat SL1, n(SL1)
*matrix SL2 = W2*Y
*svmat SL2, n(SL2)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X) (rho1:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
*ml init rho2:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo M3 
drop SL11
*/

* 2. Alternative policy (Co2 Allowances) weighted by geographical distance w/o focal policy (Co2 tax) weighted by geographical distance

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1  sigma
tempvar A rSL1  
gen double `rSL1'=`rho1'*SL1
*gen double `rSL2'=`rho2'*SL2
scalar p1 = `rho1'
*scalar p2 = `rho2'
matrix p1W1 = p1*W1
*matrix p2W2 = p2*W2
matrix IpW = I_n - p1W1
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`mu', 0, `sigma'))
end

**********************
*Open Data For Weights
**********************

* pick  the following spatial weights
spatwmat using "/Users/.../gkm_replication_isq/carbon policy/distance_gleditch_ind.dta", name(W1)  standardize

**************************
*Open Data for Regression
**************************
drop _all
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear

global nobs = 187  
matrix I_n = I($nobs)
global Y co2taxnormal
global X lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration dpi ratio_over_eutob allowances_gdp_lag
global Z allowances_gdp_lag

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W1*Z
svmat SL2, n(SL2)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y=$X SL2) (rho1:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
*ml init rho2:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max, difficult

eststo M3b


* 3. Focal policy (CO2 tax) weighted by geographical distance + Alternative policy (Co2 Allowances) weighted by geographical distance X globalization

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1  sigma
tempvar A rSL1 rSL2 
gen double `rSL1'=`rho1'*SL1
*gen double `rSL2'=`rho2'*SL2
*gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
*scalar p2 = `rho2'
*scalar p3 = `rho3'
matrix p1W1 = p1*W1
*matrix p2W2 = p2*W2
*matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1   
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`mu', 0, `sigma'))
end

**********************

* eigenval(E)
spatwmat using "/Users/.../gkm_replication_isq/carbon policy/distance_gleditch_ind.dta", name(W1)  standardize

**************************
*Open Data for Regression
**************************
drop _all
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear

global nobs = 187 
matrix I_n = I($nobs)
global Y co2taxnormal
global X lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration  dpi ratio_over_eutob allowances_gdp_lag
global Z allowances_gdp_lag

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)

matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W1*Z
svmat SL2, n(SL2)

mkmat SL2, matrix(G)

* Note: SL2 stands for (Wgeographic distance*CO2 Allowances)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************

ml model lf splag_ll (mu: $Y= lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita  co2ktpercapitasq energydependencektofoileq govteffectiveness euintegration dpi allowances_gdp_lag  c.SL2##c.ratio_over_eutob ) (rho1:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
*ml init rho2:_cons=0
*ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max

eststo M4 

margins, at(SL2 =(0(250)2000) ratio_over_eutob =(-.45)) post
estimates store marginplot1

************
* Figure 5 * 
************

clear
pr drop _all
set more off
set matsize 800

***********************
*Likelihood Evaluator  
***********************                                                         

program define splag_ll

args lnf mu rho1  sigma
tempvar A rSL1 rSL2 
gen double `rSL1'=`rho1'*SL1
*gen double `rSL2'=`rho2'*SL2
*gen double `rSL3'=`rho3'*SL3
scalar p1 = `rho1'
*scalar p2 = `rho2'
*scalar p3 = `rho3'
matrix p1W1 = p1*W1
*matrix p2W2 = p2*W2
*matrix p3W3 = p3*W3
matrix IpW = I_n - p1W1   
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`mu', 0, `sigma'))
end

**********************

* eigenval(E)
spatwmat using "/Users/.../gkm_replication_isq/carbon policy/distance_gleditch_ind.dta", name(W1)  standardize

**************************
*Open Data for Regression
**************************
drop _all
use "/Users/.../gkm_replication_isq/carbon policy/EU neighbors_carbon policy.dta", clear

global nobs = 187 
*258 if capopen is not in
matrix I_n = I($nobs)
global Y co2taxnormal
global X lag_co2taxnormal gdppercapita gdppercapitasq  co2ktpercapita  co2ktpercapitasq energydependencektofoileq govteffectiveness   euintegration  dpi    ratio_over_eutob allowances_gdp_lag
global Z allowances_gdp_lag

mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)

matrix SL1 = W1*Y
svmat SL1, n(SL1)
matrix SL2 = W1*Z
svmat SL2, n(SL2)

mkmat SL2, matrix(G)

************************
*Produce starting values
************************ 

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

***************************
*Estimate spatial lag model
***************************


ml model lf splag_ll (mu: $Y= lag_co2taxnormal gdppercapita gdppercapitasq co2ktpercapita  co2ktpercapitasq energydependencektofoileq  govteffectiveness euintegration dpi allowances_gdp_lag  c.SL2##c.ratio_over_eutob) (rho1:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
*ml init rho2:_cons=0
*ml init rho3:_cons=0
ml init sigma:_cons=`OLSsigma'
*ml trace on
ml max

margins, at(SL2 =(0(250)2000) ratio_over_eutob =(.45)) post
estimates store marginplot2

cd "/Users/.../gkm_replication_isq/carbon policy/Figures"

graph set window fontface "Times New Roman"

coefplot  (marginplot1,  lpattern(ldash) lcolor(gs13))  , at ytitle(Marginal Effect of Economic Flows) xtitle(Spatially Weighted CO2 Trade Allowances) title("Low Levels of Economic Flows") recast(line) levels(90) lwidth(*2) ciopts(recast(rline) lpattern(dash)) ylabel(-25(25)60,nogrid) xsca(alt) yline(0, lcolor(red)) plotregion(lcolor(black)) scheme(s2mono) graphregion(fcolor(white))   saving(marginplot1, replace)
coefplot  (marginplot2,  lpattern(ldash) lcolor(gs13))  , at ytitle(Marginal Effect of Economic Flows) xtitle(Spatially Weighted CO2 Trade Allowances) title("High Levels of Economic Flows") recast(line) levels(90) lwidth(*2) ciopts(recast(rline) lpattern(dash)) ylabel(-25(25)60,nogrid) xsca(alt) yline(0, lcolor(red)) plotregion(lcolor(black)) scheme(s2mono) graphregion(fcolor(white)) saving(marginplot2, replace) 

twoway histogram  SL2 if ratio_over_eutob<-0.0099, freq gap(5) color(grey14) ylabel(0(30)60,nogrid) plotregion(lcolor(black)) xlabel(0(500)2000,nogrid) fysize(22)  xtitle(Spatially Weighted CO2 Trade Allowances)  graphregion(fcolor(white)) saving(hist1, replace) 
twoway histogram  SL2 if ratio_over_eutob>-0.0099, freq gap(5) color(grey14)  ylabel(0(30)60,nogrid) plotregion(lcolor(black))  xlabel(0(500)2000,nogrid) fysize(22) xtitle(Spatially Weighted CO2 Trade Allowances)  graphregion(fcolor(white))  saving(hist2, replace) 

graph combine marginplot1.gph hist1.gph, hole(2 4) imargin(0 0 1 0) graphregion(margin(l=10 r=10)) graphregion(fcolor(white))  saving(pic1, replace)
graph combine marginplot2.gph hist2.gph, hole(2 4) imargin(0 0 1 0) graphregion(margin(l=10 r=10)) graphregion(fcolor(white))  saving(pic2, replace)  

gr combine pic1.gph pic2.gph, ycommon col(2) plotregion(fcolor(white)) graphregion(fcolor(white))   title(`"Dependent Variable: CO2 Tax Levels"', color(black))  saving(Figure5, replace)


